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ABSTRACT 

Exploring the bulge region of our Galaxy is an interesting but challenging quest because of its complex structure and the highly 
variable extinction. We re-analyse photometric near-infrared data in order to investigate why it is so hard to reach a consensus on the 
shape and density law of the bulge, as witnessed in the literature. The apparent orientation of the bulge seems to vary with the range 
of longitude, latitude, and the population considered. To solve the problem we have used the Besancon galaxy model to provide a 
scheme for parameter fitting of the structural characteristics of the bulge region. The fitting process allows the shape of the bulge's 
main structure to be determined . 

We explore various parameters and shapes for the bulge population, based on Ferrer's ellipsoids, and fit the shape of the inner disc in 
the same process. The results show that the main structure has a standard triaxial boxy shape with an orientation of about 13° with 
respect to the Sun-Galactic centre direction. But the fit is greatly improved when we add a second structure, which is a longer and 
thicker ellipsoid. We emphasize that our first ellipsoid represents the main boxy bar of the Galaxy and that the thick bulge population 
could be either (i) a classical bulge slightly flattened by the effect of the bar's potential, or (ii) an inner thick disc counterpart. 
With Ferrer's ellipsoid, the model shows a general agreement with 2MASS data at the level of 10% in the whole bulge region but does 
not produce the "double clump" feature. However, we show that the double clump seen at intermediate latitudes can be reproduced 
by adding a slight flare to the bar. To characterize the populations better, we further simulate several fields that have been surveyed in 
spectroscopy and for which a metallicity distribution function (MDF) are available. The model agrees well with these MDF measured 
along the minor axis if we assume that the main bar has a mean solar metallicity and the second thicker population has a lower 
metallicity. It then naturally creates a vertical metallicity gradient by mixing the two populations. 

In the process of model fitting, we also determine the thin disc parameters. The thin disc is found to have a scale length of 2.2 kpc, 
in good agreement with previous estimates towards the anticentre, but with a large hole of scale length 1.3 kpc, giving a maximum 
density in the plane for this population at about 2.3 kpc from the Galactic centre. In the very central part of the bulge on top of our two 
populations, and by subtracting the fitted ellipsoids, we find evidence of an extra population in the nuclear region, at 2° in longitude 
and 1° of latitude from the Galactic centre. Its location corresponds well to the central molecular zone, and to the Alard nuclear bar. 

Key words. Galaxy: stellar content. Galaxy: evolution. Galaxy: bulge. Galaxy: formation. Galaxy: dynamics 



1. Introduction 

Since the discovery of a triaxial structure in the central regions 
of th e Galaxy from 24 micron observation s (Blitz & SpergelL 
[199 1 ), COBE NIR data ('Weiland et al.',T994';'Dwek etalll995h 
or visible photometry (Staneket al., 1994), numerous attempts 
have been made to characterize this structure and to investigate 
its origin. It is still unclear whether this structure, often called the 
outer bulge but sometimes the bar, had its origin in the early for- 
mation of the spheroid (like a typical bulge, similar to ellipsoidal 
galaxies) or was formed by a bar instability later in the disc. It 
is crucial to investigate this formation history as a benchmark 
in understanding formation of disc galaxies. Kinematical studies 
are very important for tracing the stellar motions and deducing 
a dynamical history of the populations in the inner regions. New 
surveys in the near future will help solve this question, particu- 
larly the Gaia mission of the European Space Agency. 

In the meantime photometric surveys, especially in the 
infrared, can help to determine the shape and density laws 
of the stellar populations in place. The question of the 
shape, orientation and extent is still open, because one can 



find for the angle of the major axis direction with re- 
gard to the Sun-Galactic centre direction values between 10° 
and 45°, depending on the fields considered, the wavelength 
and the method use^. Among al l jhe attem£ts, the signif- 
icant ones in clude [Sta nek et al j (Il994l). il3\yek et al ( 1995|, 
IFuxI (Il997h. iLopez ^orredoira et all ^20001). iFreudenreichI 
(ll998h.lPicaud & RobinI (|2004. iBabusiaux & Gilmord (l2005h . 
Be njamin et~lj ([2005), andRattenburv et alj (l2007h . The diffi- 
culty in accurately measuring the orientation of the bulge might 
be aggravated if, actually, there is more than one population 
prese nt in the region and if these population s have different ori- 
gins. [Martinez- Valpuesta & GerhardI (1201 Ih also offer a clever 
explanation for the lack of consensus concerning the bar angle, 
which they think is partlyfrom the cone effect and partly from 
the presence of leadi ng or traiUng overdensi t ies at the end of the 
bar, as also shown bv lRomero-Gomez et alJ (1201 ih . 

In this context, metallicity distribution studies can be useful 
con straints fo r distinguishing populations. Zoccali et al. (200^ 
and lHill et all (201 1) highlight a bimodal metallicity distribution 
in Baade's window and several fields on the bulge minor axis. 
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thus identifying a mixture of populations. Concerning the kine- 
matics, proper motions are difficuh to obtain due to the large 
distance of the bulge, while radial velocities ar e expensive to 
process for large numbers of stars. ISumi et aP (l2004l) studied 
proper motions from the OGLEIII survey and compared them 
with a dynamical model. However, the OGLE fields avoid the 
Galactic plane and the northern part of the bulge, they are lim- 
ited by extinction(because the observations were performed in 
the visible), and the mean proper motion error is too large to 
accurately determine velocities at the bulge distance. 

iRich et all (l2007h have undertaken a wide survey (BRAVA, 
for bulge radial velocity assay) of radial velocities at longi- 
tudes between -10° and 10° and at latitudes -4° and -8°. Their 
first results from the latitude -4° data set show that the slope 
of the rotation curve flattens considerably at longitude \l\ > 
5° indicating a probable change in the dominant population 
at this point. Furthermore, the analy sis of metallicity versus 
kinematics by iBabusi aux et all (1201 Q t) for RGB stars in three 
bulge fields, at different latitudes along the minor axis, show 
that two populations with distinct chemo-dynamical signatures 
seem to be present in these fields. Since 2010, there have been 
several studie s showing double clumps in a number o f fields 
in the bulge ('Nat afet ali I20T0I: iMcWilUam & Zoccalil I20T0I: 
[Saito et al., 2011). This is also a feature that should be taken 
into account in an overall view of the bulge region and that the 
models should explain. 

Complete understanding of the bulge's structure and forma- 
tion would require studying the abundances and kinematics in 
various fields at different longitudes and latitudes, such as the 
APOGEE project (Majewski et al, 2007). In the mean time we 
propose to explore the bulge structure based on an analysis of 
existing photometric data. We intend to check whether a detailed 
modelling of the stellar populations and interstellar extinction 
present in the central region might lead to an answer We use the 
Two Micron All Sky Survey, hereafter 2MASS ( Skrutski e"et al.L 
I2OO 6) photometry and a model-fitti ng procedure b ased on the 
Besangon Galaxy model (BGM) dRobin et al.Ll2003l) and explore 
a wide region -20° < Z < 20° and -10°< ^ < 10° to fit a multi- 
dimensional space of parameters defining the density laws of the 
bulge and disc populations. The analysis proceeds in several 
steps. First we attempt to fit the whole area with a Monte Carlo 
scheme to explore the parameter space using the differential star 
counts in bins of K and J-K in a set of selected windows. We 
compare the resulting model with integrated star counts in order 
to visually check the realism of the resulting model. In a second 
step we compare simulated CMDs in a sample of directions to 
check their agreement in more detail. At this stage we introduce 
a modulation of the density distribution of the bar structure to 
reproduce the double clump structure seen in medium latitude 
dir ections bv .Nataf et al. (2010), McWilliam & Zoccah (2010) , 
and lSaito et al.l (l20Ilh . In a final step we compare model simula- 
tions of the metallicity distribution functions with data over the 
minor axis of the bar, placing constraints on the mean metallicity 
of the composite populations in the bulge region. 

In section 2 we describe the basis of the BGM and the free 
parameters which we adjust to the existing data. The data are 
described in section 3 along with the method that corrects for the 
3D extinction distribution towards the bulge. The fitting method 
and results are given in section 4. Validation of our extinction 
method is provided in section 5. In section 6, we compare model 
predictions with CMDs in several directions and explain how the 
model can reproduce the double clump feature in the CMDs. In 
sect. 7, we zoom in the internal bulge. In sect. 8 we discuss the 



results, and the metallicity distribution function, and present the 
perspectives of this work. 

2. Population synthesis model for the bulge 

The stellar population model is based on a population synthe- 
sis scheme. Four distinct populations are assumed (a thin disc, 
a thick disc, a bulge, and a spheroid), each deserving specific 
treatment. In the central regions, only the bulge and thin disc 
are important. The thick disc is a minor component, and the halo 
completely negligible at the magnitudes and latitudes considered 
here. The thick disc contributes slightly to the star density but 
only when the bulge and the thin disc become less important, 
that is, at latitudes above about 8-10°. 

2.1. The thin disc population 

The thin disc contributes significantly to star counts in the 
Galactic central region, so its characteristics have to be fitted 
at the same time as the bulge parameters. The important param- 
eters are the scale length and the size of the central hole. Other 
parameters and structures, such as outer radius, war p, and flare, 
are taken into account in the BGM and described in lRevle et al.l 
(2009), but are not relevant here. 

A standard evolution model is used to produce the disc pop- 
ulation, based on a set of evolutionary tracks, a constant star for- 
mation rate (hereafter SFR) over 10 Gyr, and a two-slope initial 
mass function 0(M) -Ax -M'" with a=\.6 for M< IM0 and 
Q'=3.0 forM>lM0. The preliminary tuning of the disc evolution 
pa rameters again s t releva nt observational data was described 
in 'Haywood et aJ ( Il997allbl) and later changes are explained in 
Robin etal. (2003). 

With the evolution model we populate the thin disc divid- 
ing it into seven age components. For each subcomponent, the 
distribution in absolute magnitude and effective temperature is 
obtained, assuming a star formation history constant between 
and 10 Gyr The star counts are computed using the standard 
equation of stellar statistics from the distribution in My and a 
spatial density distribution. It is therefor equivalent to assuming 
that the SFR is roughly the same over all the thin disc. 

Th e thin disc density distribution model follows the lEinastol 
(Il979h law: the distribution of each disc component (except for 
the very young one of age less than 150 million years, which 
is not important here) is described by an axisymmetric ellipsoid 
with an axis ratio depending on the age. The density law of the 
ellipsoid is described by the subtraction of two functions: 



Pd = Prfo X [exp(- ^0.25 + (— )2) - exp(- ^0.25 + (— )2)] 

with - + [^f, where: 

- R and Z are the cylindrical Galactocentric coordinates; 

- e is the axis ratio of the ellipsoid. Values of e as a function o f 
age and local normalization are given in lRobin et alJ (l2003h : 

- R^i is the scale length of the disc and is around 2.2-2.5 kpc 
(Ruphy et al. 1996). This parameter will be fitted in the pro- 
cedure described in this paper; 

- Rh is the scale length of the hole. It is also fitted here. The 
maximum density of the disc population is approximately at 
2 kpc from the centre; 

- The normalization pd„ is deduced from the local luminos- 
ity function (JahreiB et al, private communication), assum- 
ing that the Sun is located at R0=8 kpc and Z0=15 pc. 
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Alternative values for the sun-Galactocentre distance Rq are 
also considered. 

2.2. The outer bulge 

iPicaud & RobinI (12004 undertook a detailed analysis of the 
outer bulge stellar density and luminosity function by fitting 
model parameters to a set of 94 windows in the outer bulge 
situated at -S° < I < 10° and -4° < b < 4°. The data w ere 
obtained by the DENIS survey team ('Ep chtein et all Il997h us- 
ing Ks magnitude and J-Ks colour distributions. Using a Monte 
Carlo method to explore an 1 1 -dimensional space of bulge and 
disc density model parameters and a maximum likelihood test, 
they showed that the bulge follows a boxy exponential or a boxy 
sech^ profile like the one fitted to the DIRBE integrated flux in 
the ne ar infrared (JF reudenreich, 1998). These are Ferrers ellip- 
soids (lFerrerslll877h . Their resulting triaxial bulge has a major 
axis pointing towards the first quadrant with an angle of about 
10° with respect to the Sun-Galactic centre direction. A full de- 
scription of the parameter values of the bulge density law can be 
found in Picaud & Robin (2004). 

He re we reconsider the density law fitted bv lPicaud & RobinI 
(12004 . The Freudenreich S function is: 

ps = po sech^(-/?s) 

with 

7;fii^[|l|C. + |I|C,]C,/c,^|Z|C„ 
xq ya zo 

(X,Y, Z) are the cartesian coordinates in the referential of the tri- 
axial structure (X being the major axis, Y the second axis and the 
Z the third). The parameters cy and are important for explor- 
ing a wide range of shapes, from "disky" to "boxy". This allows 
great flexibility: one can even have a "disky" shape in the plane, 
together with a boxy projection vertically. 

In the following we also consider alternative functions for 
the overall shape, exponential (E), and Gaussian (G): 

Pe = Po exp(-7?i) 

and 

PC = po exp(-;?j)^ 

The density function is then multiplied by the cut-off func- 
tion fc (distances given in kpc, and Rc is called the cut-off" ra- 
dius): 

P = P fciRxY), with RxY = VX2 + F2 
RxY ^ Rc fc(RxY) = 1 

RxY >Rc^ MRxy) = exp (-(^)') 
Three angles define the orientation: 

- (p: orientation angle from the sun-centre direction of the pro- 
jection on the Galactic plane of the bulge major axis, 

- j6: pitch angle of the bulge major axis from the Galactic 
plane, 

- y: roll angle around the bulge major axis. 

In the iPicaud & RobinI (l2004 study, the angle /3 was found 
very close to 0°, and the third angle y was found to be ill defined 
because the minor axes Y and Z had similar scale lengths. In this 
study, we therefore adopt j6=y=0°. Finally, the angle 0, the three 



scale lengths xq, yo, zo, the density at the centre po, the cut-off" 
radius R^ and the two coefficients C\\ and Cj_ are considered as 
free density parameters in the fitting process. 

The bulge stellar population is taken from a single burst 
population of ages var ying between 6 and 10 Gyr as tested by 
IPicaud & Robin[(l2004 . The favoured combinations between the 
age and evolution models are from Bruzual et a D([i997) with an 
age of 10 Gyr, or from Girardi et al (2002) with an age of 7.94 
Gyr (log(age)=9.9). It appears that the Kj band counts are not 
very sensitive to the age, so a proper estimation of the age and 
age range of the populations in the bulge would require comple- 
ment ary data. Thi s analysis thus takes the luminosity function 
from lGirardi et all (12002) with a log(age) of 9.9 as a reference. 

2.3. One or two populations in the bulge? 

In the bulge region there may be several populations, with sev- 
eral star formation histories, cohabiting. The coexistence of dif- 
ferent populations in the bulge region can be tested using the 
distribution in colour-magnitude diagrams. In our fitting process 
we consider alternatively one or two bulge populations, and for 
each of them the S, E, or G functions. 

3. Data description and extinction corrections 

3.1. 2MASSclata 

The IPicaud & RobinI (|2004 analysis was limited in longitude 
and latitude (-12°< I < -4° < b <4°), used DENIS data, 
and did not contain many fields near the plane (less than 10% 
were situated at \l\ < 1°). Here, we use the 2MASS point source 
catalogue data. We concent rate on the region -20° < I <20°, 
-10° < b <10°. Unlike the iPicaud & RobinI (|2004 study, the 
selection goes down to latitude b=0°. The data are used in three 
ways: (i) for d etermining the extinction using the method from 
Marshall et all (,2006 ). (ii) for fitting bulge model parameters on a 
selection of windows, (iii) for comparing the whole bulge region 
with the fitted model and checking the extinction distribution 
again. 

3.2. Extinction 

Since the extinction has a major eff'ect on star counts close to 
the plane, a detailed analysis of the extinction field by field 
is needed to provide realistic model sim ulations. For this pur - 
pose, we used the method described in Marshall et all (l2006l) . 
This method adopts the standard BGM and determines the ex- 
tinction by comparing the simulated colour-magnitude diagrams 
from this model and the observational data. Adjusting both the 
extinction and the bulge parameters in a single run is a difficult 
data analysis problem. Therefore we proceeded with iterations, 
first fitting the extinctio n with the preliminary bulge model (from 
IPicaud & RobinI |2004 . then adjusting the bulge parameters as- 
suming this extinction and finally verifying that the extinction 
derived from the new fitted bulge model is not significantly dif- 
ferent from the original extinction distribution. 

The goodness of fit used in the extinction determination is 
based on a chi-squared statistic for two bin ned distributions hav- 
ing diff'erent total number of elements ( P ress et al 1 119921) . As 
such, the extinction determination is sensitive to diff"erences in 
colour and not to discrepancies between the number of observed 
and modelled stars. 

In the present study, the Ks magnitude and the J-Kj colour 
were used to compare with simulations. Cuts were made in Ks 
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Fig. 1. Distribution in longitude and latitude of the 200 windows 
used in the fitting process. The colour coding corresponds to the 
limiting magnitude of the counts in Kg . 



and J at the completeness limit, defined from the histogram in 
magnitude at the bin before the maximum for limiting the risk of 
incompleteness that can occur even at the peak in magnitude. In 
the extinction determination and in the maps for visualizing the 
comparison between model and data (see section 4. 1 and figures 
2, 3, and 4) we further restrict the comparison at Kj <12 because 
the extinction determination method is only based on the giants' 
colour and that they can be separated from the dwarfs only up to 
this magnitude. 



4. Fitting procedure 

The first set of data we use is a selection of fields from the 
2MASS point source catalogue. Two hundred windows have 
been chosen in order to cover the region defined by -20° < 
/ <20°and -10°< b <10°. For each window, the extinction is 
determined using the method described in section 3.2, as well as 
the completeness limit independently in each band. The BGM in 
its standard version is used to simulate the NIR data, taking the 
completeness into account. To illustrate, we show in fig. [T] the 
200 windows and the completeness limit in the Kj band. Then 
the following procedure is used to adjust the bulge model pa- 
rameters. 

For each window the star counts in the CMD diagram are 
binned in Kj and J-K^ in bins with equal numbers of stars, two 
bins in colour, and seven in magnitude, giving a total of 14 bins. 
In the model the counts of each population are counted sepa- 
rately to be able to change them as a function of the variable 
parameters. The parameter space is explored randomly from suc- 
cessive Monte Carlo drawings and a selection of the best models 
in the sen se for the likelihood. T he likelihood function is the 
same as in Pi caud & Robinl(l200? ). The process is stopped when 
the variation in the parameters from one iteration to the next is 
less than the required accuracy for each parameter The process 
is done 20 to 50 times to ensure a robust solution. 

In a second step, we compare the fitted model with 2MASS 
data in the whole region -20°< / < 20°and -10°< b < 10° 
and compute the relative difference between data and model on 
a grid of 15x15 arcmin fields. The aim is twofold. Firstly, it en- 
sures that the 2MASS data are correctly fitted by the model all 
over the bulge area (not only in the 200 selected windows), with- 
out systematics in any region. Secondly, it verifies that the new 
extinction, recomputed from the new density model, has not var- 



ied significantly since the density fitting process, so we ensure 
that proceeding in two separate steps (fitting of the extinction, 
then fitting of the density model) is efficient and will converge. 

We first attempted to fit a unique bulge population all over 
the considered region. The parameters that are fitted are the fol- 
lowing: the bulge angle (angle between the major axis of the 
bulge and the Sun-Galactic centre direction), the bulge scale 
lengths (considering 3 axes, xo, yo, zq), the bulge mass, the ex- 
tent in radial distance (cutoff R^), two parameters describing the 
boxiness of the bulge cy and Cj_ (see above), the disc scale length, 
and the scale length describing the size of the hole inside the 
disc. 

As we show below, the agreement was not good enough, and 
we had to consider alternative models. We attempted to fit a sec- 
ond structure, with the same kinds of parameters as for the first 
bulge structure. We also considered different formulae for the 
bulge density and for the second structure: sech^, Gaussian, or 
exponential. 

We assume here a Sun to G alactocentric distan ce of 8 
kpc (following recent results from ' Ghez et al] (1200 8^). unlike 
Picaud & Robin (2004) who assume a value of 8.5 kpc. We in- 
vestigate the effect of this value on the result in the discussion. 

In the next subsections we describe the results for each kind 
of model fitted. 

4.1. Fitting a single ellipsoid 

The fit of a single ellipsoid on the large area considered here 
leads to s ignificantly diffe r ent res ults compared with the analysis 
done by iPicaud & RobinI (l2004l) where the data set considered 
was restricted to the region -12°< / < 8° and -4°< b < 4° with 
about 100 windows. The important thing here is probably not the 
number of windows, which is doubled, but rather the extent of 
the area. Moreover, unlike the Picaud & Robin (2004) study, we 
drop the selection of giants and consider a fit covering the whole 
colour range. Finally, we account for interstellar extinction in a 
more realistic manner than they did. 

Results are given in Table 1 . Model A indicates the parame- 
ters fitted on DENIS data by Picaud & Robin (2004) for compar- 
ison only. In the model A fitting process, the parameters cy and 
c± were fixed to 2 and the Sun-galactocentric radius to 8.5 kpc. 
For each model and each parameter, the dispersion around the 
fitted values is indicated in the second line. When a model has 
a parameter fixed rather than fitted, it is in italics. The B model 
is the model with a free angle with the (default) sech^ shape. C 
is similar but with a fixed angle of 20°, D is similar but with a 
fixed angle of 25°. The G model is similar to the B model but 
with a Gaussian shape, the E model has an exponential shape 
(see below). 

The most significant difference between the new result 
(model B) and the old one (model A) is the new scale length 
of the bulg e (major axis scale lengt h xq), which was found to be 
1.59 kpc in'Picaud & RobinI (12004 ) and is here about 4 kpc. The 
cutoff radius has also grown from 2.54 to about 6 kpc. It means 
that the limited region used in our first study has led to under- 
estimating of the bulge's extent. The second axis scale length 
is also greater than in Picaud & Robin (2004) (corresponding to 
the depth of the bulge along the line of sight if the angle is small), 
while the scale height (minor axis Zq) has not changed (the lati- 
tude range covered by the previous paper was sufficient to cover 
the bulge). We also found that the scale length of the disc has 
slightly diminished and that the hole has diminished. 

We plot in figure |2] star counts in Kj from 2MASS data (top 
left) compared with the fitted model with single population (top 
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Table 1. Best fit models with 1 single bulge ellipsoid. The columns give the fitted parameters and reduced likelihood Lr (fi: bulge 
angle of the major axis with the Sun-Centre distance in degrees; xO, yO, zO: scale lengths along the 3 axes in kpc; normalization 
factor; Rc: cutoff radius in kpc, cy and : parameters of boxyness, Rd: disc scale length in kpc, Rh: disc hole scale length in kpc. 



Model 





xO 


yO 


zO 


Normalization 


Rc 


C|| 


c± 


Rd 


Rh 


Lr 






kpc 


kpc 


kpc 


l.e9Mopc-3 


kpc 






kpc 


kpc 




A 


11.10 


1.59 


0.424 


0.424 


9.63 


2.54 


2. 


2. 


2.53 


1.32 




B 


7.1 


4.07 


0.76 


0.41 


23.83 


5.99 


1.434 


3.797 


2.26 


0.18 


-22842. 


stdev 


1.01 


0.263 


0.043 


0.018 


1.31 


0.025 


0.162 


0.858 


0.165 


0.190 


892. 


C 


20.0 


4.48 


1.46 


0.47 


26.47 


5.98 


1.015 


4.579 


2.64 


0.99 


-32895. 


stdev 





0.403 


0.122 


0.023 


0.928 


0.0188 


0.135 


0.491 


0.302 


0.183 


1707. 


D 


25.0 


3.77 


1.48 


0.46 


29.09 


6.00 


1.019 


4.379 


2.47 


1.03 


-40027. 


stdev 





0.325 


0.084 


0.024 


1.09 


0.022 


0.123 


0.542 


0.321 


0.283 


1578. 


G 


5.8 


4.12 


0.84 


0.41 


16.43 


5.99 


1.271 


2.223 


2.30 


0.61 


-22943. 


stdev 


1.12 


0.260 


0.121 


0.016 


0.273 


0.039 


0.197 


0.870 


0.208 


0.250 


609. 


E 


6.9 


2.54 


0.45 


0.23 


49.75 


5.98 


1.476 


4.656 


2.25 


0.60 


-23661. 


stdev 


1.21 


0.310 


0.050 


0.013 


2.85 


0.083 


0.381 


0.977 


0.167 


0.247 


1514. 



right), residuals as defined by (A^mod - A^obs)/A^obs (bottom left), 
and difference divided by the Poissonian counting error of the 
model (i.e. (Nobs-N,nod)/ V(^mod) (bottom right). When comput- 
ing the extinction, the Kj band was limited to 12 or brighter. The 
other bands have no such limitation. This limit ensures lower 
uncertainty on the 2MASS PSC fluxes. This figure shows the 
B model described in Table 1 after the 3D extinction distribu- 
tion has been refitted. The structure of the residuals shows that 
the shape of the fitted bulge model is not appropriate. In the in- 
plane region the model density agrees well with the data. But the 
fields at high latitudes and intermediate longitudes have X-shape 
residuals, where the model density is too low. It is probable that 
the resulting fit is a bad compromise between adjusting the den- 
sities in plane and at different latitudes. It is worth noting that in 
the top left-hand panel the 2MASS star count map indeed shows 
boxy iso-densities that are not reproduced by the model, even 
if the density laws were chosen expUcitly to allow for a boxy 
shape. 

Interestingly, the residual map shows a structure in the inner 
bulge at longitudes less th an 2°, which we believe is the nuclear 
bar seen bV iAlardI (1200 Ih . This population has not yet been in- 
cluded in the model, and is shown here by subtraction. We dis- 
cuss the point in section 7. 

4.2. Influence of Other parameters 

The Sun-Galactocentric distance Rq is not adjusted here dur- 
ing the fitting process. It is fixed at either 7., 7.5, 8, or 8.5 kpc. 
Models with Rq of 8 kpc give the larger likelihood, but the data 
are not sensitive enough to distances along the line of sight to 
strongly constrain this parameter. Moreover, changing Rq has 
little impact on the other parameters, such as cy and Cj^ , which 
stay of the order of 1 and 4 (resp). 

We assume an age of about 8 Gyr for the bulge population 
log(ag e)=9.9), following the best result in the iPicaud & RobinI 
20041) study. Our analysis is not very sensitive to the assumed 
age for an old population (age over 5 Gyr) in the red clump. 
The difference in the luminosity function would be noticeable 
if we had covered the red clump and the giant branch down to 
the turnoff of the population, which is not possible with 2MASS 



data. The turnoff position and the giant sequence are also in- 
fluenced by the age. The turnoff would be slightly brighter in 
a younger population (but this is not visible in most 2MASS 
fields), and the giant sequence would be slightly bluer But this 
effect is of second order and diflicult to see in broad band pho- 
tometry and even more in extincted fields. Complementary data 
would be necessary to constrain the age of the population at 
the level of a Gyr Consequently, we consi der an isochrone fo r 
the bulge with an age of about 8 Gyr from iGirardi et all (l2002h . 
which should be valid if the true age of the population is in 
the range 6-10 Gyr Alternative isochrones within a reasonable 
range of age would not significantly change the conclusion on 
the bulge shape and the extinction. 

4.3. About the bulge/bar angle 

Many bulge studies find a significantly larger bulge angle than 
our result (from one population model B : 7.1°). Most of them 
are in the range 10° to 45°. We found 11.1° in model A. Since 
some of our parameters are slightly correlated, we attempted to 
fix the bulge angle at different values to see whether it is possible 
to find a good enough fit with a larger bulge angle and to see in 
which field this value is favoured. Models C and D in Table 1 
have fixed bulges angle of 20° and 25° respectively. These mod- 
els have very low likelihood with regards to smaller angles. 

We suspect that the bulge angle found in different stud- 
ies might depend on the field selection, extent in longi- 
tude, distance to the Galactic plane, and depth of the survey. 
Mar tinez- Valpuesta & GerhardI (1201 Ih have shown how the bar 
angle determination can be strongly biased by the cone volume 
effect, in the case of the long bar, but this effect is also present 
when determining the main bulge angle using the mean distance 
of the red clump. 

To test these effects, we selected subsets of fields in longitude 
and latitude and reproduced the fitting procedure. It appears that 
for fields at |/| > 12° and close to the plane (\b\ < 3°), the likeli- 
hood is higher with a large angle, while this large angle gives a 
bad fit at low longitude and high latitude (\b\ > 3° and |Z| < 5°). 
Fields at high longitudes and low latitudes are well fitted by large 
bulge angles, while fields at low longitudes and high latitudes 
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Fig. 2. Star counts up to limiting magnitude Ks from 2MASS data (top left) compared with fitted 1 ellipsoid model B (top right), 
residuals (Nmod-Nobs)/ Nobs (bottom left), and difference divided by the Poissonian counting error of the model (i.e. (Nobs-Nmod)/ 
V(Mnod) (bottom right). In the residual map, contours are drawn at intervals of 20% model overestimate (black) and 20% model 
underestimate (white). The 1 -ellipsoid model leaves significant X-shaped residuals. Near the Galactic centre the nuclear bar popu- 
lation is missing in the model. The residuals in the outer regions are not very significant due to the small number of stars in each bin 
(seen in the top left panel in dark violet and black). 



favour small bulge angles. When fixing the angle of the bar at 
higher values (models C and D, 20 or 25°), the fields situated 
at large longitudes show similar likelihood, but fields at higher 
latitudes (\b\ > 3°) get signi ficantly worse. This is i n good agree- 
ment with the analysis from lCabrera-Lavers et alj (|2007, 2008), 
from 2MASS data as well, but limited to fewer fields. These find- 
ings might indicate that a single ellipsoid fit is not adequate. This 
is why in the following section we attempt to fit the sum of two 
ellipsoids on our data set. 

4.4. Fitting two ellipsoids 

We attempt to apply the procedure to fit two ellipsoids (in addi- 
tion to the thin disc) having the same luminosity function, age, 
and metallicity, as these parameters do not influence the over- 
all near infrared photometry very much (see the discussion be- 
low). The same parameter set describing each ellipsoid is fit- 
ted on both populations. We consider several combinations of 
sh apes, S shap e as before (S), exponential (E), Gaussian (G) as 
in lDweket air(ll995l) . 

Table 2 shows the parameters of the best fit with different 
shapes S, E, or G. Because it is unlikely that two populations 
have different angles from a dynamical point of view, and be- 
cause the fit with diff'erent angles did not give significantly dif- 
ferent angles, we have imposed the two ellipsoids to have the 
same orientation in the plane. 

Among the best combinations obtained for the best model 
shapes, the parameters are very similar: the best fit is obtained 
with the two ellipsoids having an orientation of about 9° to 13°, 
the first ellipsoid with a scale length Xq of 1.6 and secondary 
scale length of 0.5 and 0.4 kpc. Values are slightly different 
when an exponential shape is assumed, because the exponential 
is significantly more like a disc than a sech^ or a Gaussian, and 
the different shape is compensated for by different scale lengths 
and normalization. Regarding the likelihood, the best solution 
is with a sech^ for the first component and an exponential for the 
second one. 



The second ellipsoid is generally less massive than the first 
one but with much larger scale lengths (4./1.3/0.8). The disc 
scale length is about 2.2 to 2.3 kpc with a wide hole of about 
1 to 1.2 kpc (but 0.5 kpc only when an E-shape is assumed for 
the bulge). As a result, the maximum of the disc density appears 
to be outside the bulge/bar maximum, as if the internal disc was 
swept out by the bar. The boxiness of the two ellipsoids is no- 
ticeable (c|| and generally larger than 2) but their values are 
not very well constrained. 

Maps of the star counts in the K, band are shown in fig- 
ure [3] to be compared with fig. |2l Residuals are much smaller 
than in the 1 -ellipsoid model, and the fit appears good in nearly 
all regions with residuals at the level of 10%. The population 
of the nuclear bar at |/| < 2° as in previous model appears by 
subtraction. To see where the differences are most significant, 
accounting for the Poisson noise, we have plotted the differ- 
ence divided by the Poissonian counting error of the model (i.e. 
(A^obs - A^mod) / V-^mod) in the bottom right-hand panel of figure[3] 

4.5. Attempts with other shapes 

As is the case for many disc galaxies, the Milky Way has a spi- 
ral structure even if it is not well defined at present. The spiral 
arms are mainly detected in the visible from young populations 
and from interstellar matter tracers such as HII regions. 2MASS 
star counts are dominated by K giants, which generally do not 
follow the spiral structure because of their age, but the counts 
could be slightly sensitive to them because in any case, a part 
of these K giants are indeed young objects. Another structure 
that has been well studied in CO is the molecular ring where 
a large proportion of the giant molecular clouds reside. It is not 
well established whethe r the stellar populations contribute to this 
ring. lSevenster & Kalna is (2001) propose that a stellar ring sig- 
nificantly enhances the microlensing optical depth towards the 
bulge. 

To check this point we investigate the model likelihood when 
a ring shape is included in the fit in place of the second ellipsoid. 
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Table 2. Parameters and reduced likelihood Lr when two ellipsoids are fitted, including main sequence and giants, and forcing the 
two ellipsoids at the same orientation. Parameters are the same as in previous table. The last two lines of the table give the standard 
deviation about the mean for model (S+E). Standard deviation for other models are similar 
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Fig. 3. Same as figure|2]but with two fitted ellipsoids. 

For this we use the ISevenster & KalnaisI (1200 ll) elliptical ring 
density law, with free parameters: ring mean radius, ellipticity, 
and orientation in the plane, thickness in the plane, scale height, 
and normalization. We also tested cutoffs in longitude for this 
ring to allow for the existence of a portion of ring rather than 
a complete one. The result is that the likelihood is comparable 
to the 1 ellipsoid model, and the density of such a structure is 
found to be compatible with 0. It means that in the longitude and 




20 10 -10 -20 

I [deg] 




20 10 -10 -20 

I [deg] 



latitude zones considered here, the ring and the spiral arms do 
not significantly contribute to the 2MASS star counts. 

4.6. Uncertainties and correlations 

Using our method we are able to estimate both the accuracy 
with which the parameters are determined and the correlations 
between them. A Monte-Carlo process was computed about 40 
times in the case of the best model (at least 20 times in other 
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Fig. 4. Distribution of 
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solutions for the main bulge population (first ellipsoid) with regards to its parameters. Columns are, from the 
in degrees, scale length xO, yO, zO in kpc, central mass, cy and . Rows are given in the same order from top 
gns indicate all likelihood, squares hold for the 10 best likelihoods, the blue circle for the best solution. 



cases) to compute the average best parameters and give estimates 
of the uncertainties and correlations between parameters. For the 
model S+E, correlations between parameters are shown in fig- 
ure |4] for the first ellipsoid and in figure |5] for the second one. 
Cross-correlations between parameters of the two ellipsoids are 
shown in figure |6] In these figures for all trials, the best ten so- 
lutions (in the sense of the likelihood) and the best likelihood 
solution shown in the table. The correlations between parame- 
ters of a given population are rather weak in general, indicating 
that the method is robust. The only net correlation appears be- 
tween scale lengths for the first ellipsoid (between xq and yo and 
Zo) and with the normalization. 

The uncertainties on the scale lengths is about 250 pc for the 
major axis, and only 50 pc on the secondary axis. The normal- 
ization is determined at the level of 15% for the main ellipsoid 
and the second ellipsoid. The boxyness parameters cy and Cj_ are 
nearly always found to be larger than two, being three in the 
mean, making both ellipsoids rather boxy, even if the secondary 
one is found with an exponential shape (which is more like a disc 
than the sech^). 

The disc and hole-scale lengths are slightly anti-correlated, 
a longer scale length favouring a smaller hole. This is expected 
due to the uncertainty on the distance indicators in broad band 
photometry. A better fit of the disc scale length would need to in- 
clude larger longitudes. But then the spiral structure should also 
be included, as we would reach arm tangents. The model also 



has a slight tendancy to produce too many stars at Z < -5° and 
too few stars at Z > 5° close to the plane. We could interpret this 
by the ability of the ho le to be ellipsoidal in the plane, as found 
by i Binnev et al.l (1 19971) . We attempted to model this ellipticity 
but it is poorly constrained, and the fit is only very slightly better 
than with a round hole. This case could be considered later in 
more detail with deeper data. 



5. Extinction validation 

The analysis was done assuming a 3D extinc tion distribution 
model obtained using the bulge model from iPicaud & RobinI 
(l200l . We need to ensure that the extinction resulting from 
the new fit is consistent with the initial extinction we assumed. 
Figure|7]shows the comparison between the projected extinction 
map used originally with the map deduced from the fit with the 
new model. The differences are negligible. We also compared the 
distribution in distance of Ak along some lines of sight and saw 
generally good agreement in the distances of the clouds, within 
the error bars of this distance (a few hundred parsecs). Figure |8] 
shows the variation in the distribution of extinction along four 
lines of sight for the standard model (old), the model with the 
bulge population simulated by one ellipsoid (Fl), and the model 
with two ellipsoids (model 2, F2). The diff'erences are generally 
negligible. In one case (longitude 0), we see that the cloud is 
found at a greater distance with the new model, but the total ex- 
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Fig. 5. Distribution of solutions for the second ellipsoid with regards to its parameters. Red plus signs indicate all solutions, squares 
hold for the 10 best likelihood solutions, blue circle for the best solution. 
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Fig. 6. Distribution of solutions and correlation between parameters of the second ellipsoid with regard to the first ellipsoid. The first 
ellipsoid is in abscissa, second in ordinate. Red plus signs indicate all likelihoods, green squares hold for the 10 best likelihoods, 
blue circle for the best solution. 



tinction in the bulge remains the same, so the extinction deduced 
using our method ( Marshall et al, 2006) is sufficiently robust and 
not particularly model dependent. Conversely, we can emphasize 
that the 3D extinction model is reliable enough to avoid biasing 
our results concerning the stellar populations study. 



6. Resulting colour-magnitude diagrams and 
double-clump feature 

We present colour-magnitude diagrams for five directions in 
order to see the subtle diff'erences between 1 -ellipsoid and 2- 
ellipsoid models. To reproduce the natural width of the se- 
quences in the CMD, we have applied a star-by-star variation 
of the extinction, about the mean extinction at the star distance 
given by the extinction map. The cosmic variation in the extinc- 
tion is assumed to be 20% of the extinction, a value that is esti- 
mated by trial and error. 

Figuresl9l fT0lTn fT2l andOshow these CMDs. The 2MASS 
data are in the first columns, 1 -ellipsoid model in column 2, 2- 
ellipsoid model in columns 3, and histograms in Ks in column 5. 
In column 4 we also present a modified model that produces a 
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Fig. 7. Ak map deduced from Marshall et al. (2006) method. 
Top: based on the standard stellar population model from BGM; 
Bottom: based on the new stellar population model. 
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Fig. 8. Distribution of extinction over four lines of sight (lati- 
tudes b=-3° and longitudes between -3 and 0°) from the orignal 
model (Old), 1 -ellipsoid model (Fl), 2-ellipsoid model (F2). 




Fig. 14. Luminosity functions for the red clump region in fields 
at latitude -8°. Coordinates in longitude and latitude are indi- 
cated in each panel. The double-clump feature is simulated by 
adding a slight flare to the original bar model (see text). 



double-clump feature (see below). One notices the incomplete- 
ness of the data starting probably at about K=12.5 in fields at 
1^1 < -4°. 

The 1 -ellipsoid model gives star counts in the red giant 
branch that are too high. The agreement is generally much bet- 
ter with the 2-ellipsoid model, even if it is not perfect in some 
fields. There are still residuals in the shape of the giant branch 
that could probably be solved by varying the star by star dis- 
persion in the extinction, especially in the case of the field with 
coordinates (l,b) - (5; 4.5). The important point here is to have 
a general model reproducing most fields, at the expense of some 
discrepancy in a limited number of fields, since we are interested 
in the large-scale structure. 

In fields at intermediate latitudes {\b\ > 5°) several stud- 
ies have shown there are double clumps. These double clump s 
change as a function of longitude (iMcWilUam & ZoccaliLl20Tob . 
From our analysis limited to 2MASS dififerential star counts, and 
due to the lack of completeness of the data, our solution does 
not present double clumps, but just a modulation of the shape 
of the luminosity function due to the shape of the bar. The sec- 
ond structure has its clump peak at fainter magnitudes than the 
2MASS completeness limit. If there is a structure creating the 
double clump it should be in the main bar. 

Since our model reproduces the counts at the level of 10% 
in the whole bulge region, we can consider that what is creating 
the double clump feature is a modulation of the whole shape 
of the bar Pfenniger ( 1985) did a numerical study of instability 
in the 3D family of periodic orbits in a barred galaxy model. He 
shows a family of resonance creating the growth of a box-shaped 
bulge. It shows how stars are diff'used at high z. It is natural to 
hypothesize that the double clump can be produced by such an 
effect. 

To test this hypothesis, we simulated a flaring bar by enhanc- 
ing the scale height of the bar with a sinus fonction of the dis- 
tance along the major axis of the bar. If x is the galactocentric 
distance along the major axis, and dzD is the scale height along 
the z axis, the new vertical scale height is simulated by 

dz^dzOxil+ 0.3 X sin(x/650.). 

The original value of dzD (obtained by the fit above) was de- 
creased slightly (from 390 to 330 pc) to account for the fact that 



we measured the mean scale height when we fixed this param- 
eter in the fit. The parameters in the equation are a reasonable 
guess in order to fit the clump data at a medium latitude. 

In figure[T4]we show the star counts around the clump in var- 
ious field s at latitu de b=-8° in order to compare with figure 3 of 
Mc Willi am & Zocc ali (2010). It shows that a flaring bar as mod- 
elled here can qualitatively reproduce the double-clump feature. 
The proportion of stars in each clump varies as a function of 
the longitude due to the flare, and the significant differences be- 
tween positive and negative longitudes. This figure can be com- 
pared with figure 3 of McWilliam & Zoccali (2010). It can also 
be seen on a zoom of the luminosity function in the field at 1=0, 
b=-7°(figure [TSl l where we superimposed the observed counts 
with the model predictions without the flare and with the flare. 
We show that a simple modulation of our primary model with a 
slight bar flare reproduces the two-clump feature. 

Since the 2MASS data are not complete, we cannot perform 
a quantitative fit of the flare parameters with these data. We only 
show here that a flaring bar can be a realistic explanation for 
the double-clump an d that it is dynamically reasonable, as in 
the iPfennige^ (1 19851) simulations. In a future paper we plan a 
more complete analysis of the shape of the flare using VVV data, 
for exam ple, by comparing simulations with the density plots 
from Sai to et al.l(l201lb in order to specify the characteristics of 
this flare more accurately as a function of the galacto-centric 
distance. 

7. The inner bulge 

Assuming that the model of stellar populations and 3D extinc- 
tion now reproduces the populations of the outer bulge and the 
disc, we can investigate if any further structures are missing from 
our modelling. As a check we plot in figure[T6]the higher reso- 
lution map of the difference between model predictions and ob- 
served star counts in the region -4°<1<4° and -3°<b<3°. The 
excess of stars in the data should represent the po pulation of the 
inner bulge. It resembles the population seen bv lAlardl(l200lh 
and coincide with the position of the central molecular zone. 
Several authors (Morris & Serabvn, 1996; An et al., 2011) have 
shown that this region contains young stellar populations. To in- 
vestigate further this population, one would need to improve the 
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Fig. 9. Colour-magnitude diagrams for the 2MASS field at 1=0, b=-4. (a) data; (b) best-fit model with 1 ellipsoid; (c) best-fit models 
with 2 ellipsoids; (d) modified model with flared bar; (e) histograms of data (red solid) and models (1 ellipsoid: green long dashed; 
2 ellipsoids; blue dotted, flared bar: magenta short dashed). 








Fig. 10. Colour-magnitude diagrams for the 2MASS field at l=-i-3, b=-3. Same coding as in fig|9l 






0.5 1 1.5 
J-K 





Fig. 11. Colour-magnitude diagrams for the 2MASS field at 1=5, b=4.5. Same coding as in fig|9l 



3D extinction model here. As we saw in figure |7] the extinc- 
tion determined by our method saturates at about A]i-3 due to 
the limiting magnitude in 2MASS. Thus our extinction is most 
probably underestimated, implying that the density of the inner 
bulge population is overestimated. We do not investigate further 
the density of this population in this paper, as it could only be 
quaUtative from 2MASS data alone. We envisage to apply our 
method to deeper data and at larger wavelengths, like Spitzer at 
3.6 and 4.5 microns, in order to make a deeper 3D map of the 
extinction in this region and to study the inner bulge region, the 



probable nuclear bar or ring and its hnk with the central molec- 
ular zone. 



8. Discussion 

We have shown in this analysis that the bulge region can be sim- 
ulated using the sum of two ellipsoids: a main component that is 
a quite standard boxy bulge/bar, which dominates the counts up 
to latitudes of about 5°, and a second ellipsoid that is a thicker 
structure seen mainly at higher latitudes. This thick component 
seems to be about as thick as the local thick disc. It could be 
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Fig. 12. Colour-magnitude diagrams for the 2MASS field at 1=-17, b=-4.25. Same coding as in fig|9l 








Fig. 13. Colour-magnitude diagrams for the 2MASS field at 1=0, b=-8. Same coding as in figH) In this field we see clearly in the 
histogram the double clump feature reproduced by the modified model (panel d and magenta line) and not by the originally fitted 
model (panel c and blue line). 



Fig, 15. Luminosity functions for the red clump region in a field 
at latitude -7°and longitude 0° from 2MASS (red dots with error 
bars) compared with the model without flare (in green solid line) 
and with the flaring bar (in blue dotted line). 
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0.50 

0.00 

-0.50 
-1.00 



Fig. 16. Difference map for the inner bulge region (Nobs- 
Nmod)/Nobs- The population in excess in the data is most likely 
the nuclear bar. 



either the classical bulge, flattened by the effect of the barred 
potential, or the inner counter part of the thick disc population. 
To distinguish how these structures can be differentiated and to 
understand which population they belong to, it is useful to look 
at the spectroscopic surveys showing metalhcity distributions at 
different latitudes. 



8.1. Metallicity as a function of latitude 

Among studies measuring the metallicity distribution function 
(MDF) in different fields of the bulge region several have con- 
sidered the characteristic s of the populations along the minor 
axis at different latitudes. IZoccah et al] (l2008h st udied the MDF 
and radial velocity distribution from the RGB, iGonzalez et al.l 
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(1201 ih completed the data with the abundances of alpha abun- 
dances, iBabusiaux et al I (l2010t) studied the correlation between 
metallicity and kinematics and FHill et al (2011) computed the 
MDF from the red clump giants for Baade's window. The over- 
all conclusion of these studies is that the MDF looks bimodal 
and the contribution of the components changes with distance 
from the plane. At a low latitude (-4°), the dominant popula- 
tion has a metallicity slightly above-solar and a small contribu- 
tion comes from a seco nd component of metallicity around -0.5. 
IBabusiaux et al I (l2010l) show that the high-metallicity compo- 
nent has a high vertex deviation compared to the low-metallicity 
component that does not. They conclude that the vertical metal- 
licity gradient is t he effect of mixing of populations. On the other 
hand. lNess et al.l ([201 1) decompose the bulge populations from 
an analysis of the ARGOS survey in up to five components to re- 
produce the complete MDF. They find a smaller g radient ( about 
0.4 dex/kpc) compared with IZoccali et al.l (l2008h . who found a 
mean gradient of 0.6 dex/kpc. 

Can we explain these MDFs by our two-ellipsoid model? 
To a nswer this question, we simulated the observed MDF from 
IZoc cali et al. (2008) in three latitude ranges, on the minor axis 
and deduced what would be the necessary metallicity of each of 
our components to explain the observations. 

We assume that the boxy bar has a mean solar metallicity 
and that the thick component has a mean metallicity of -0.35 
dex. We did not fit the density of each population, which is the 
one obtained from the 2MASS star-count fitting. W e applied the 
selecti on functio n expla ined in lZoccali et al.l (l2008h for the RGB 
and in iHill et all ( 1201 ll) for the red clump sample. Figure [ITj 
shows the comparison between the observed metallicity distri- 
bution and the simulated one for RGB giants at l=0°and b=-4°, 
-6°, and -12°, and for red clump giants at b=-4°. At the bot- 
tom we also show the decomposition of the distribution in the 
three major components: the thin disc, the "boxy-bar", and the 
"thick-bulge". We see good agreement beween our model pre- 
dictions and the data in all four cases, and even clearer at high 
latitudes, within the small Poisson statistics of the samples. The 
error bars are only the Poisson noise on the counts in the selected 
sample. It does not take observational errors on the metallicity 
into account or the noise due to the selection function. The only 
significant disagreement comes from a lower number of low- 
metallicity stars in the model for the field at b=-6°compared to 
the data. This feature will be explored in more detail as soon as 
a larger sample will be available. It will then be possible to de- 
termine whether a metallicity gradient should be added to one 
of the two populations or whether the relative density of the two 
populations is represented well by this model or should be re- 
vised given the errors. 

8.2. Conclusion and future plans 

We have presented here an attempt to fit differential star counts, 
CMDs, and metallicity distribution functions in a wide area over 
the Milky Way bulge region. The model accounts for the metal- 
licity distribution on the minor axis of the bulge, the dissym- 
metry of the star counts as a function of longitude, and the ex- 
istence of double clumps at medium latitudes. This model has 
two components. The main one is modelled by boxy, S-shape 
triaxial Ferrer ellipsoids, having scale lengths of 1.46/0.49/0.39 
kpc, an angle of the main axis with regards to the Sun-Galactic 
centre direction of 13°, and a total mass of 6.1 x IO'^Mq . The 
age is assumed to be 8 Gyr or more, and the mean metallicity 
is approximately solar. The second component is modelled by a 
triaxial exponential Ferrer ellipsoid less massive (2.6 x 1O**)M0 , 



with scale lengths of 4.44/1.31/0.80 kpc, also old with a poorer 
metallicity of -0.35 dex. A slight flare of the bar component 
is able to explain the double clumps qualitativeley at latitudes 
higher than 6°. 

This model with two components explains the observed 
metallicity distribution function and its gradient along the mi- 
nor axis by a transition between the two components dominat- 
ing the counts, and also better reproduces the detailed distribu- 
tion in the CMD s. A preliminary comparison with Brava data 
iRich et a DdlOOT') shows that the main component is a fast rota- 
tor, so is probably a boxy bar (the object of a paper in prepara- 
tion). The second structure (the thicker bulge) has a smaller ve- 
locity dispersion, thus could be either a classical bulge flattened 
by the potential of the bar or formed by early mergers, as in the 
iBournaud et al.l (l2007h scenario. Testing of these scenarios are 
on-going, using our model and comparing predictions with kine- 
matical and metallicity data in the whole bulge region. Th e forth- 
coming surveys, APOGEE proiect lMaiewski et all (l2007h . Gaia- 
ESO spectroscopic survey, among others, and later the ES A Gaia 
mission ( http://www.rssd. esa.int/index.php?project=GAIA) are 
all expected to offer new insights into the kinematics and abun- 
dances of the central region of the Galaxy, giving strong con- 
straints on any scenario for the formation of the bulge and bar 
stellar populations. We emphasize that our model can be used 
for preparing these future surveys and will be a useful tool for 
their interpretation. 
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Bottom: decomposition in populations of the model simulation. Red solid line: total; green long dash: disc, short dash blue: boxy 
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